classdef SimpleNMRG < handle
    %SIMPLENMRG �˴���ʾ�йش����ժҪ
    %   �˴���ʾ��ϸ˵��
    
    properties
        session
        device
        zi
        
        channel129
        channel131
        demodCh129
        demodCh131
        channelMainCurrent
        main_current_control_volt
        
        frequency129
        frequency131
        amplitude129
        amplitude131
        
        resonance129    % record the latest resonance frequency and other info of Xe129, every time we get a new data, please update this var.
        resonance131    % TODO: verify this update
    end
    
    methods
        function obj = SimpleNMRG(session, device, varargin)
            %SIMPLENMRG ��������ʵ��
            %   �˴���ʾ��ϸ˵��
            L=log4m.getLogger;
            
            p=inputParser;
            isSession = @(x) isa(x, 'MeasurementSession.AbstractMeasurementSession');
            chTest = @(x) ismember(x, {'0', '1', '2', '3'});
            ctrlVoltTest = @(x) 0<=x && x <=10;
                        
            p.addRequired('session', isSession);
            p.addRequired('device', @ischar);
            p.addParameter('channel129', '0', chTest);
            p.addParameter('channel131', '1', chTest);
            p.addParameter('demodCh129', '2', chTest);
            p.addParameter('demodCh131', '3', chTest);            
            p.addParameter('amplitude129', 0.015, @isnumeric);
            p.addParameter('amplitude131', 0.040, @isnumeric);
            p.addParameter('main_current_channel', '0', chTest);
            p.addParameter('main_current_control_volt', 8.948, ctrlVoltTest);
            p.parse(session, device, varargin{:});      
            
            obj.session = session;
            obj.device = device;
            obj.zi = obj.session.getDevice(device);

            obj.channel129 = p.Results.channel129;
            obj.channel131 = p.Results.channel131;
            obj.demodCh129 = p.Results.demodCh129;
            obj.demodCh131 = p.Results.demodCh131;
            obj.amplitude129 = p.Results.amplitude129;
            obj.amplitude131 = p.Results.amplitude131;
            obj.channelMainCurrent = p.Results.main_current_channel;
            obj.main_current_control_volt = p.Results.main_current_control_volt;
            
            %
            obj.zi.loadSettings( [Constant.AutoMeasurementCorePath '\+EquipmentControl\+SystemlMachine\@SimpleNMRG\' 'settings_HighFieldNMR_Matlab.xml']);

            %demod setting
            obj.zi.setInputSignalSelect(obj.demodCh129, 0); % 0 - Signal in
            obj.zi.setDemodOscSelect(obj.demodCh129, obj.channel129);   % for signal analyse
            obj.zi.setDemodOscSelect(obj.channel129, obj.channel129);   % for driving output
            obj.zi.setDemodPhase(obj.demodCh129, 0.0);
            obj.zi.setDemodHarm(obj.demodCh129, 1);
            
            obj.zi.setInputSignalSelect(obj.demodCh131, 0); % 0 - Signal in
            obj.zi.setDemodOscSelect(obj.demodCh131, obj.channel131);   % for signal analyse
            obj.zi.setDemodOscSelect(obj.channel131, obj.channel131);   % for driving output
            obj.zi.setDemodPhase(obj.demodCh131, 0.0);
            obj.zi.setDemodHarm(obj.demodCh131, 1);
            
            obj.zi.disableOutAmplitude(); % disable all driving
            obj.zi.setAuxOutSelect('0', -1);
            obj.zi.setAuxOutSelect('1', -1);
            obj.zi.setAuxOutSelect('2', -1);
            obj.zi.setAuxOutSelect('3', -1);
            

            %default resoance
            res129.frequency = 258.100; res129.amplitude = 0.120; res129.gamma = 0.050; res129.phase = [];
            res131.frequency =  76.500; res131.amplitude = 0.020; res131.gamma = 0.050; res131.phase = [];
            obj.resonance129 = res129;
            obj.resonance131 = res131;
            
            L.info('SimpleNMRG', 'Object created.');
        end
        %%
        function initDriving(obj, varargin)
            L=log4m.getLogger;
            
            p=inputParser;
            p.addParameter('mainCurrentControlVolt', 8.948);
            p.addParameter('nmr129freq', 258.126);
            p.addParameter('nmr131freq', 76.525);
            p.addParameter('amp129', 0.015);
            p.addParameter('amp131', 0.040);
            p.addParameter('drvOffset', 0.100);
            p.parse(varargin{:});

            obj.resonance129.frequency = p.Results.nmr129freq;
            obj.resonance131.frequency = p.Results.nmr131freq;
            obj.setMainCurrent(p.Results.mainCurrentControlVolt)
            obj.setDrivingFrequency129(p.Results.nmr129freq).setDrivingAmplitude129(p.Results.amp129);
            obj.setDrivingFrequency131(p.Results.nmr131freq).setDrivingAmplitude131(p.Results.amp131);
            obj.setDrivingOffset(p.Results.drvOffset).nmrDrivingOn();
            fprintf('wainting 30s...');
            pause(30);
            fprintf('Done.\n');
            obj.setZeroPhase129().setZeroPhase131();
            pause(5);
            
            L.info('SimpleNMRG::initDriving', 'Driving condition initialized.');
        end

        %% NMR Resonance
        function data=findResonance(obj, varargin)
            expectedName={'129', '131', 'both'};
            p=inputParser;
            p.addParameter('nuclearName', '129', @(x) any(validatestring(x,expectedName)) );
            p.parse(varargin{:});

            data.fid=obj.fidSignal('nuclearName', p.Results.nuclearName);
            data.fft=obj.fftSignal(data.fid,'nuclearName', p.Results.nuclearName);
            data.fit=obj.fitSignal(data.fft,'nuclearName', p.Results.nuclearName);

        end
        
        function data=fidSignal(obj, varargin)
            expected_nuclearName={'129', '131', 'both'};
            p=inputParser;
            p.addParameter('nuclearName', '129', @(x) any(validatestring(x,expected_nuclearName)) );
            p.addParameter('sampling_rate', 12, @isnumeric);
            p.addParameter('frame_time', 60.0, @isnumeric);
            p.addParameter('addToSession', true, @islogical);            
%             p.addParameter('trig_level_129', 0.070, @isnumeric);            
%             p.addParameter('trig_level_131', 0.010, @isnumeric);           
            p.parse(varargin{:});
            
            L=log4m.getLogger;
            L.info(['SimpleNMRG::fidNMR',p.Results.nuclearName], ['NMR', p.Results.nuclearName, 'FID measurement started.']);
            clear textprogressbar;
            
            scp=obj.zi.scopeModule();
            scp.enableChannel('1');
            scp.inputSelect('0', 0);           
            scp.reset();
            scp.setSamplingRate(p.Results.sampling_rate); % index = 14 : 3.66 kHz;  12 : 14.6kHz
            scp.setFrameLength(p.Results.frame_time);
            scp.disableTrigger().setContinuousAcquisition(p.Results.frame_time);
            
%             scp.setTriggerSlope(2); % fall
%             scp.setTriggerHysteresis(0.0);
%             scp.disableTriggerGating();
%             scp.setTriggerHoldoff(0.0);
%             scp.setTriggerReference(0.0);
%             scp.setTriggerDelay(0.0);           
%             scp.enableTrigger().setTimeOut(10.0).setSingleShot();
            
            scp.subscribeChannel();
            
            switch p.Results.nuclearName
                case '129'
%                     scp.setTriggerChannel(50); % Trigger Channel 50 = Demod3R
%                     scp.setTriggerLevel(p.Results.trig_level_129);
                    obj.disableDriving131();
%                     fprintf('%s: Waiting for Pre-relaxation (%d seconds).\n', datestr(now,'HH:MM:SS.FFF, mm/dd/yy'), 120);
%                     pause(120);
                    obj.disableDriving129();                   
                    fid=scp.run();                    
                    obj.enableDriving129();
                    obj.enableDriving131();
                case '131'
%                     scp.setTriggerChannel(51); 
%                     scp.setTriggerLevel(p.Results.trig_level_131);
                    obj.disableDriving129();
                    fprintf('%s: Waiting for Pre-relaxation (%d seconds).\n', datestr(now,'HH:MM:SS.FFF, mm/dd/yy'), 120);
                    pause(120);
                    obj.disableDriving131();     
                    fid=scp.run();
                    obj.enableDriving131();
                    obj.enableDriving129();
                case 'both'                    
                    obj.disableDriving129();     
                    obj.disableDriving131();     
                    fid=scp.run();
                    obj.enableDriving129();
                    obj.enableDriving131();
            end
            
            if p.Results.addToSession
                scp.addToSession(['fidNMR', p.Results.nuclearName]);
            end
            scp.unsubscribeChannel();
            
            data=fid.getNodeData;
            
            figure('Name', 'FID')
            plot(data.time, data.value, '.b', 'linewidth', 1.5);
            grid on;
            title(['fidNMR', p.Results.nuclearName]);
            xlabel('Time (s)');
            ylabel('Amplitude');
            ylim([-1.1*max(abs(data.value)), 1.1*max(abs(data.value))]);

            L.info(['SimpleNMRG::fidNMR', p.Results.nuclearName], ['NMR', p.Results.nuclearName, 'FID measurement finished.']); 
        end
        
        function data=fftSignal(obj, fidData, varargin)
            expected_signalName={'Re', 'Im', 'Abs'};
            expected_nuclearName={'129', '131', 'both'};
            p=inputParser;
            p.addParameter('seleSignal', 'Abs', @(x) any(validatestring(x, expected_signalName)) );
            p.addParameter('nuclearName', '129', @(x) any(validatestring(x, expected_nuclearName)) );
            p.addParameter('bandwidth', 1.0, @isnumeric);
            p.addParameter('estFrequency129', 258.2, @isnumeric);
            p.addParameter('estFrequency131', 76.5, @isnumeric);
            p.parse(varargin{:});
            
            samplingRate=1/(fidData.time(2)-fidData.time(1));
            fidData.value=double([fidData.value;zeros(2^24-length(fidData.value),1)]);
            points=length(fidData.value);
            data.frequency=samplingRate*(0:points/2)/points;
            switch p.Results.seleSignal
                case 'Re'
                    data.amplitude=2*real(fft(fidData.value))/points;
                case 'Im'
                    data.amplitude=2*imag(fft(fidData.value))/points;
                case 'Abs'
                    data.amplitude=2*abs(fft(fidData.value))/points;
            end           
            data.amplitude=data.amplitude(1:points/2+1);
            
            switch p.Results.nuclearName
                case '129'
                    resonance_range=find(data.frequency>=p.Results.estFrequency129-p.Results.bandwidth & data.frequency<=p.Results.estFrequency129+p.Results.bandwidth);
                case '131'
                    resonance_range=find(data.frequency>=p.Results.estFrequency131-p.Results.bandwidth & data.frequency<=p.Results.estFrequency131+p.Results.bandwidth);
                case 'both'
                    bandwidth2=p.Results.estFrequency129-p.Results.estFrequency131;
                    resonance_range=find(data.frequency>=p.Results.estFrequency131- 0.1*bandwidth2 & data.frequency<=p.Results.estFrequency129+0.1*bandwidth2);
            end
            data.frequency=data.frequency(resonance_range);          
            data.amplitude=data.amplitude(resonance_range)/max(data.amplitude(resonance_range));
            [vMax, idx]=max(abs(data.amplitude));
            
            figure('Name', ['fftNMR',p.Results.nuclearName])
            plot(data.frequency, data.amplitude, '.b', 'linewidth', 1.5);
            grid on;
            title(sprintf('maxFind: w_0^{129}=%3.6f Hz', data.frequency(idx)));
            switch p.Results.nuclearName
                case 'both'
                    [~, idx129] =max(data.amplitude);
                    [~, idx131] =max(data.amplitude(1:end/2));
                    title(sprintf('maxFind:  w0=%3.6f Hz for NMR131; w0=%3.6f Hz for NMR129', data.frequency(idx131), data.frequency(idx129)));
                otherwise
                    title(sprintf('maxFind: w0=%3.6f Hz for NMR%s', data.frequency(idx), p.Results.nuclearName));
            end
            xlabel('Frequency (Hz)');
            ylabel('Normalized Amplitude');
            ylim([0, 1.1*vMax]);
            
        end
        
        function data=fitSignal(obj, fftData, varargin)
            expected_nuclearName={'129', '131', 'both'};
            p=inputParser;
            p.addParameter('nuclearName', '129', @(x) any(validatestring(x,expected_nuclearName)) );
            p.addParameter('bandwidth', 1.0, @isnumeric);
            p.addParameter('estFrequency129', 258.2, @isnumeric);
            p.addParameter('estFrequency131', 76.5, @isnumeric);
            p.parse(varargin{:});
            
            [xData, yData] = prepareCurveData( fftData.frequency, fftData.amplitude );
            [vMax, idx] =max(abs(yData));           
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            
            switch p.Results.nuclearName
                case 'both'
                    ft=obj.fitModel('modelName', 'SqrtLorentz2');
                    [~, idx129] =max(yData);
                    [~, idx131] =max(yData(1:end/2));
                    opts.Lower =      [0.0  0.0  0.0  0.0    xData(idx131)-p.Results.bandwidth, 0.0, xData(idx129)-p.Results.bandwidth]; % a, b, c, gamma131, w0_131, gamma129, w0_129
                    opts.Upper =      [Inf  Inf  Inf  1.0    xData(idx131)+p.Results.bandwidth, 1.0, xData(idx129)+p.Results.bandwidth];
                    opts.StartPoint = [1    1    0.0  0.01   xData(idx131), 0.01, xData(idx129)];
                otherwise
                    ft=obj.fitModel('modelName', 'SqrtLorentz');                    
                    opts.Lower =      [0.0  0.0  0     min(fftData.frequency)]; % a, b, gamma, w0
                    opts.Upper =      [Inf Inf 1.0   max(fftData.frequency)];
                    opts.StartPoint = [vMax 0.0  0.010 fftData.frequency(idx)];
            end

            [data.fitResult, data.gof] = fit( xData, yData, ft, opts );
            
            data.frequency=fftData.frequency;
            data.amplitude=feval(data.fitResult, data.frequency);
            
            figure('Name', ['fitNMR',p.Results.nuclearName])
            plot(data.fitResult,fftData.frequency,fftData.amplitude);
            grid on;
            xlabel('Frequency (Hz)');
            ylabel('Normalized Amplitude');
            ylim([0, 1.1*vMax]);
            switch p.Results.nuclearName
                case '129'
                    title(sprintf('fitcoeffs: w0=%3.6f Hz, tau=%2.1f s for NMR%s', data.fitResult.w0, 1/(2*pi*data.fitResult.gamma), p.Results.nuclearName));
                    obj.resonance129.frequency = data.fitResult.w0;
                    obj.go2Resonance129(20);
                case '131'
                    title(sprintf('fitcoeffs: w0=%3.6f Hz, tau=%2.1f s for NMR%s', data.fitResult.w0, 1/(2*pi*data.fitResult.gamma), p.Results.nuclearName));
                    obj.resonance131.frequency = data.fitResult.w0;
                    obj.go2Resonance131(20);
                case 'both'
                    title(sprintf('fitcoeffs:  w0=%3.6f Hz, tau=%2.1f s for NMR131; w0=%3.6f Hz, tau=%2.1f s for NMR129', data.fitResult.w0_131, 1/(2*pi*data.fitResult.gamma131), data.fitResult.w0_129, 1/(2*pi*data.fitResult.gamma129)));                    
            end

        end
        
        function ft=fitModel(obj, varargin)
            expected_modelName={'Lorentz','Gaussian','PseudoVoigt','SqrtLorentz','SqrtLorentz2'};
            p=inputParser;
            p.addParameter('modelName', 'SqrtLorentz', @(x) any(validatestring(x, expected_modelName)));
            p.parse(varargin{:});
            
            switch p.Results.modelName
                case 'Lorentz'
                    ft=fittype('a*gamma/(pi*((w-w0)^2+gamma^2))+b', 'indep', 'w');
                case 'Gaussian'
                    ft=fittype('a*exp(-(w-w0)^2/(2*sigma^2))/(sqrt(2*pi)*sigma)+b', 'indep', 'w');
                case 'PseudoVoigt'
                    ft=fittype('a*( eta*gamma/(pi*((w-w0)^2+gamma^2))+(1-eta)*exp(-(w-w0)^2/(2*sigma^2))/(sqrt(2*pi)*sigma) )+b', 'indep', 'w');
                case 'SqrtLorentz'
                    ft=fittype('a*gamma/sqrt(pi*((w-w0)^2+gamma^2))+b', 'indep', 'w');
                case 'SqrtLorentz2'
                    ft=fittype('a*gamma131/sqrt(pi*((w-w0_131)^2+gamma131^2))+b*gamma129/sqrt(pi*((w-w0_129)^2+gamma129^2))+c', 'indep', 'w');
            end
            
        end
  
        
        function nmr129 = sweepNMR129(obj, varargin)
            L=log4m.getLogger;
            L.info('SimpleNMRG::sweepNMR129', 'NMR129 sweep started.');
            
            clear textprogressbar;

            p=inputParser; df = 0.020; %20mHz
            p.addParameter('start', obj.resonance129.frequency-df, @isnumeric);
            p.addParameter('stop', obj.resonance129.frequency+df, @isnumeric);
            p.addParameter('samplecount', 50, @isnumeric);
            p.addParameter('bandwidth', 22.0, @isnumeric);
            p.addParameter('order', 8, @isnumeric);
            p.addParameter('settling_time', 1, @isnumeric);
            p.addParameter('settling_inaccuracy', 0.1e-3, @isnumeric);
            p.addParameter('averaging_sample', 100, @isnumeric);
            p.addParameter('averaging_time', 0, @isnumeric);
            p.addParameter('addToSession', true, @islogical);
            p.parse(varargin{:});  
            
            % set frequency to start frequency for pre-relaxation
            obj.setDrivingFrequency129(p.Results.start);
            
            swpr=obj.zi.sweepModule;
            
            swpr.reset();
            swpr.setControl('gridnode',  ['oscs/' obj.channel129 '/freq'], ...
                            'start', p.Results.start, 'stop', p.Results.stop, 'samplecount', p.Results.samplecount);
            swpr.setFilter('bandwidthcontrol', 1, ... %bandwidth contrl - fix
                           'bandwidth', p.Results.bandwidth, 'order', p.Results.order);
            swpr.setSettling('settling_time', p.Results.settling_time, 'settling_inaccuracy', p.Results.settling_inaccuracy);
            swpr.setStatistics('averaging_sample', p.Results.averaging_sample, 'averaging_time', p.Results.averaging_time);

            swpr.setMonitorScale('linear', 'linear');
            swpr.subscribeChannel(obj.demodCh129, 'r');
            
            % wait for pre-relaxation
            fprintf('%s: Waiting for Pre-relaxation (%d seconds).\n', datestr(now,'HH:MM:SS.FFF, mm/dd/yy'), 20);
            pause(20);
            
            nmr129 = swpr.run();
            swpr.unsubscribeChannel();
            
            [freq, amp, gamma] = obj.resonance_fit(obj.demodCh129, nmr129);
            obj.resonance129.frequency = freq;
            obj.resonance129.amplitude = amp;
            obj.resonance129.gamma = gamma;          
            
            if p.Results.addToSession
                swpr.addToSession('NMR129');
            end
            L.info('SimpleNMRG::sweepNMR129', 'NMR129 sweep finished.');
        end
        function nmr131 = sweepNMR131(obj, varargin)
            L=log4m.getLogger;
            L.info('SimpleNMRG::sweepNMR131', 'NMR131 sweep started.');

            clear textprogressbar;

            p=inputParser; df = 0.200; %20mHz
            p.addParameter('start', obj.resonance131.frequency-df, @isnumeric);
            p.addParameter('stop', obj.resonance131.frequency+df, @isnumeric);
            p.addParameter('samplecount', 50, @isnumeric);
            p.addParameter('bandwidth', 22.0, @isnumeric);
            p.addParameter('order', 8, @isnumeric);
            p.addParameter('settling_time', 1, @isnumeric);
            p.addParameter('settling_inaccuracy', 0.1e-3, @isnumeric);
            p.addParameter('averaging_sample', 100, @isnumeric);
            p.addParameter('averaging_time', 0, @isnumeric);
            p.addParameter('addToSession', true, @islogical);
            p.parse(varargin{:});  
            
            % set frequency to start frequency for pre-relaxation
            obj.setDrivingFrequency131(p.Results.start);
            
            swpr=obj.zi.sweepModule;

            swpr.reset();
            swpr.setControl('gridnode',  ['oscs/' obj.channel131 '/freq'], ...
                            'start', p.Results.start, 'stop', p.Results.stop, 'samplecount', p.Results.samplecount);
            swpr.setFilter('bandwidthcontrol', 1, ... %bandwidth contrl - fix
                           'bandwidth', p.Results.bandwidth, 'order', p.Results.order);
            swpr.setSettling('settling_time', p.Results.settling_time, 'settling_inaccuracy', p.Results.settling_inaccuracy);
            swpr.setStatistics('averaging_sample', p.Results.averaging_sample, 'averaging_time', p.Results.averaging_time);

            swpr.setMonitorScale('linear', 'linear');
            swpr.subscribeChannel(obj.demodCh131, 'r');
            
            % wait for pre-relaxation
            fprintf('%s: Waiting for Pre-relaxation (%d seconds).\n', datestr(now,'HH:MM:SS.FFF, mm/dd/yy'), 20);
            pause(20);
            
            nmr131 = swpr.run();
            swpr.unsubscribeChannel();
            
            [freq, amp, gamma] = obj.resonance_fit(obj.demodCh131, nmr131);
            obj.resonance131.frequency = freq;
            obj.resonance131.amplitude = amp;
            obj.resonance131.gamma = gamma;
                        
            if p.Results.addToSession
                swpr.addToSession('NMR131');
            end
            
            L.info('SimpleNMRG::sweepNMR131', 'NMR131 sweep finished.');
        end        
        
        function obj = go2Resonance(obj, pause_time)
            L=log4m.getLogger;
            
            if nargin == 1
                pause_time = 10;
            end
            
            obj.setDrivingFrequency129(obj.resonance129.frequency);            
            obj.setDrivingFrequency131(obj.resonance131.frequency);
            pause(pause_time);
            
            obj.zi.DemodPhaseZero(obj.demodCh129);
            obj.zi.DemodPhaseZero(obj.demodCh131);
            
            obj.resonance129.phase = obj.zi.getDemodPhase(obj.demodCh129);
            obj.resonance131.phase = obj.zi.getDemodPhase(obj.demodCh131);
            L.info('SimpleNMRG::go2Resonance', 'Rosonant points of Xe129 and Xe131 reset.');
        end
        function obj = go2Resonance129(obj, pause_time)
            L=log4m.getLogger;
            
            if nargin == 1
                pause_time = 10;
            end
            obj.setDrivingFrequency129(obj.resonance129.frequency);            
            pause(pause_time);
            
            obj.zi.DemodPhaseZero(obj.demodCh129);
            obj.resonance129.phase = obj.zi.getDemodPhase(obj.demodCh129);
            
            L.info('SimpleNMRG::go2Resonance129', 'Rosonant point of Xe129 reset.');
        end
        function obj = go2Resonance131(obj, pause_time)
            L=log4m.getLogger;
            
            if nargin == 1
                pause_time = 10;
            end
            obj.setDrivingFrequency131(obj.resonance131.frequency);
            pause(pause_time);

            obj.zi.DemodPhaseZero(obj.demodCh131);
            obj.resonance131.phase = obj.zi.getDemodPhase(obj.demodCh131);
            
            L.info('SimpleNMRG::go2Resonance131', 'Rosonant point of Xe131 reset.');
        end
        
        function last_sweep_resonance(obj)
            obj.resonance_fit(obj.demodCh129, obj.session.measurement_data.NMR129(end) );
            obj.resonance_fit(obj.demodCh131, obj.session.measurement_data.NMR131(end) );
        end


        %% Main Current
        function obj = setMainCurrent(obj, volt)
            L=log4m.getLogger;
            
            if nargin == 1
                volt = obj.main_current_control_volt;
            end
            
            v0 = obj.zi.getAuxOutOffset(obj.channelMainCurrent);
            dv = min( abs(volt-v0), volt/1000.0);
            dv = sign(volt-v0)*dv;
            
            for v = v0:dv:volt
                obj.zi.setAuxOutOffset(obj.channelMainCurrent, v);
                pause(0.01);
            end
            
            L.info('SimpleNMRG::setMainCurrent', sprintf('Main current control set to %5.3f V.', volt) );
        end
        
        %% driving field
        function obj = nmrDrivingOn(obj)
            obj.zi.fireSigOut();
        end
        function obj = nmrDrivingOff(obj)
            obj.zi.silenceSigOut();
        end
        
        function obj = enableDriving129(obj)
            obj.zi.enableOutAmplitude(obj.channel129);
        end
        function obj = disableDriving129(obj)
            obj.zi.disableOutAmplitude(obj.channel129);
        end
        function obj = enableDriving131(obj)
            obj.zi.enableOutAmplitude(obj.channel131);
        end
        function obj = disableDriving131(obj)
            obj.zi.disableOutAmplitude(obj.channel131);
        end
        function obj = setDrivingAmplitude129(obj, amplitude)
            if nargin == 1
                amplitude = obj.amplitude129;
            end
            obj.zi.setSigOutAmplitude(obj.channel129, amplitude);
        end
        function obj = setDrivingAmplitude131(obj, amplitude)
            if nargin == 1
                amplitude = obj.amplitude131;
            end
            obj.zi.setSigOutAmplitude(obj.channel131, amplitude);
        end
        function obj = setDrivingOffset(obj, offset)
            obj.zi.setSigOutOffset(offset);
        end
        
        function obj = setDrivingFrequency129(obj, frequency)
            obj.zi.setOscFreq(obj.channel129, frequency);
        end
        function obj = setDrivingFrequency131(obj, frequency)
            obj.zi.setOscFreq(obj.channel131, frequency);
        end
        
        function obj = setZeroPhase129(obj)
            obj.zi.DemodPhaseZero(obj.demodCh129);
        end
        function obj = setZeroPhase131(obj)
            obj.zi.DemodPhaseZero(obj.demodCh131);
        end

        
        %% OpenLoop Data
        function data = pollData(obj, storagename, varargin)
            clear textprogressbar;
            
            p=inputParser;
            durationTest = @(x) x > 0;
            p.addParameter('duration', 5, durationTest)
            p.parse(varargin{:});
            
            plotter=obj.zi.plotterModule(obj.session);
            
            plotter.reset();
            plotter.setMonitorScale('linear', 'linear');
            plotter.subscribeChannel(obj.demodCh129, 'x');
            plotter.subscribeChannel(obj.demodCh129, 'y');
            plotter.subscribeChannel(obj.demodCh129, 'r');
            plotter.subscribeChannel(obj.demodCh129, 'theta');
            plotter.subscribeChannel(obj.demodCh129, 'frequency');
            
            plotter.subscribeChannel(obj.demodCh131, 'x');
            plotter.subscribeChannel(obj.demodCh131, 'y');
            plotter.subscribeChannel(obj.demodCh131, 'r');
            plotter.subscribeChannel(obj.demodCh131, 'theta');
            plotter.subscribeChannel(obj.demodCh131, 'frequency');
            
            data = plotter.run( 'total_duration', p.Results.duration, ...
                                'window_time', p.Results.duration);
            plotter.unsubscribeChannel();
            
            try
                plotter.addToSession(storagename);
            catch
                warning('Data not add to session.');
            end

        end
        
        %% spectrum OpenLoop
        function data = spectrumData(obj, storagename, varargin)
            clear textprogressbar;
            
            p=inputParser;
            nonNeg = @(x) x > 0;
            p.addParameter('maxF', 200, nonNeg);
            p.addParameter('nPoint', 8192, nonNeg);
            p.addParameter('nAverage', 100, nonNeg);
            p.parse(varargin{:});
            
            spect = obj.zi.spectrumModule();

            spect.setControl('grid_repetitions', p.Results.nAverage);
            spect.setParam('spectrum/frequencyspan', p.Results.maxF);
            spect.setParam('grid/cols', p.Results.nPoint);

            spect.reset();
            spect.setMonitorScale('log', 'log');
            
            spect.subscribeChannel('2', 'sample.X.fft.abs.avg');
            spect.subscribeChannel('2', 'sample.Y.fft.abs.avg');
            spect.subscribeChannel('2', 'sample.R.fft.abs.avg');
            spect.subscribeChannel('2', 'sample.Theta.fft.abs.avg');
            spect.subscribeChannel('2', 'sample.Frequency.fft.abs.avg');
            
            spect.subscribeChannel('3', 'sample.X.fft.abs.avg');
            spect.subscribeChannel('3', 'sample.Y.fft.abs.avg');
            spect.subscribeChannel('3', 'sample.R.fft.abs.avg');
            spect.subscribeChannel('3', 'sample.Theta.fft.abs.avg');
            spect.subscribeChannel('3', 'sample.Frequency.fft.abs.avg');
            
            data= spect.run();
            spect.unsubscribeChannel();

            try
                spect.addToSession(storagename);
            catch
                warning('Data not add to session.');
            end
        end
        
        %% Lock
        function obj = lock129(obj, varargin)
            L=log4m.getLogger;
            
            p=inputParser;
            p.addParameter('P', -0.002, @isnumeric);
            p.addParameter('I', -0.002, @isnumeric);
            p.addParameter('D', 0.0, @isnumeric);
            p.parse(varargin{:});
            
            ziDAQ('setInt',    ['/' obj.device '/pids/0/mode'], 1);
            ziDAQ('setInt',    ['/' obj.device '/pids/0/inputchannel'], 2);
            ziDAQ('setInt',    ['/' obj.device '/pids/0/outputchannel'], 0);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/setpoint'], 0.0);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/demod/timeconstant'], ziBW2TC(22, 8) );
            ziDAQ('setInt',    ['/' obj.device '/pids/0/demod/order'], 8);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/demod/harmonic'], 1);
            ziDAQ('setInt',    ['/' obj.device '/pids/0/phaseunwrap'], 0);
            
            ziDAQ('setDouble', ['/' obj.device '/pids/0/center'], obj.resonance129.frequency);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/limitlower'], -0.020);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/limitupper'], 0.020);
            
            ziDAQ('setDouble', ['/' obj.device '/pids/0/p'], p.Results.P);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/i'], p.Results.I);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/d'], p.Results.D);
            ziDAQ('setDouble', ['/' obj.device '/pids/0/dlimittimeconstant'], 0);

            ziDAQ('setDouble', ['/' obj.device '/pids/0/rate'], 2142000);
            
            ziDAQ('setInt',    ['/' obj.device '/pids/0/enable'], 1);
            
            L.info('SimpleNMRG::lock129', sprintf('Xe129 resonance locked: P=%f, I=%f, D=%f', p.Results.P, p.Results.I, p.Results.D) );
        end
        function obj = unLock129(obj)
            L=log4m.getLogger;
            ziDAQ('setInt',    ['/' obj.device '/pids/0/enable'], 0);
            L.info('SimpleNMRG::lock129', 'Xe129 resonance unlocked.');
        end
        
        function obj = lock131(obj, varargin)
            L=log4m.getLogger;
            p=inputParser;
            p.addParameter('P', -0.002, @isnumeric);
            p.addParameter('I', -0.002, @isnumeric);
            p.addParameter('D', 0.0, @isnumeric);
            p.parse(varargin{:});
            
            ziDAQ('setInt',    ['/' obj.device '/pids/1/mode'], 1);
            ziDAQ('setInt',    ['/' obj.device '/pids/1/inputchannel'], 3);
            ziDAQ('setInt',    ['/' obj.device '/pids/1/outputchannel'], 1);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/setpoint'], 0.0);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/demod/timeconstant'], ziBW2TC(22, 8) );
            ziDAQ('setInt',    ['/' obj.device '/pids/1/demod/order'], 8);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/demod/harmonic'], 1);
            ziDAQ('setInt',    ['/' obj.device '/pids/1/phaseunwrap'], 0);
            
            ziDAQ('setDouble', ['/' obj.device '/pids/1/center'], obj.resonance131.frequency);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/limitlower'], -0.020);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/limitupper'], 0.020);
            
            ziDAQ('setDouble', ['/' obj.device '/pids/1/p'], p.Results.P);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/i'], p.Results.I);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/d'], p.Results.D);
            ziDAQ('setDouble', ['/' obj.device '/pids/1/dlimittimeconstant'], 0);

            ziDAQ('setDouble', ['/' obj.device '/pids/1/rate'], 2142000);
            
            ziDAQ('setInt',    ['/' obj.device '/pids/1/enable'], 1);
            L.info('SimpleNMRG::lock131', sprintf('Xe131 resonance locked: P=%f, I=%f, D=%f', p.Results.P, p.Results.I, p.Results.D) );
        end
        function obj = unLock131(obj)
            L=log4m.getLogger;
            ziDAQ('setInt',    ['/' obj.device '/pids/1/enable'], 0);
            L.info('SimpleNMRG::lock131', 'Xe131 resonance unlocked.');
        end
        
        %% scope
        function data = signalScope(obj, storagename)
            clear textprogressbar;
            
            scp=obj.zi.scopeModule();
            
            scp.enableChannel('both');
            scp.inputSelect('0', 0);
            scp.inputSelect('1', 12);

            %%
            scp.reset();
            scp.disableTrigger();
            scp.setMode('frequency');
            scp.setFFTWindow('hann');
            scp.setSamplingRate(9);
            scp.setLength(65536);
            scp.setAveraging(100);
            scp.setContinuousAcquisition(60.0);
            
            scp.setMonitorScale('log', 'log');

            scp.subscribeChannel();
            data=scp.run();
            scp.unsubscribeChannel();
            try
                scp.addToSession(storagename);
            catch
                warning('Data not add to session.');
            end
        end

    end
    methods(Static)
        function [frequency, amplitude, gamma, base, fitresult, gof] = resonance_fit(channel, nmr_sweep_data)
            data = nmr_sweep_data.getNodeData(channel, 'r');
            phase_data = nmr_sweep_data.getNodeData(channel, 'phase');
            xList = data.grid; yList = data.value;
            [vMax, idx] = max(data.value);
            % Fit: 
            [xData, yData] = prepareCurveData( xList, yList );

            % Set up fittype and options.
            ft = fittype( 'A*g*g/((x-x0)*(x-x0) + g*g)+B', 'independent', 'x', 'dependent', 'y' );
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            opts.Lower =      [0.0  0.0  0     min(data.grid)]; % A B g x0
            opts.Upper =      [10.0 10.0 1.0   max(data.grid)];
            opts.StartPoint = [vMax 0.0  0.050 xList(idx)];
            

            % Fit model to data.
            [fitresult, gof] = fit( xData, yData, ft, opts );
            
            frequency = fitresult.x0;
            amplitude = fitresult.A + fitresult.B;
            gamma = fitresult.g;
            base = fitresult.B;
            
            idx = find(phase_data.grid<frequency);
            idx1 = idx(end);
            idx2 = idx1+1;
            diffFreq=phase_data.grid(idx2) - phase_data.grid(idx1);
            diffPhase=phase_data.value(idx2) - phase_data.value(idx1);
            slope=diffPhase/diffFreq/1000/pi*180;
            
            % Plot fit with data.
            figure( 'Name', 'NMR Resonance');
            subplot(2, 1, 1);
            plot( fitresult, xData, yData); grid on; hold on;
            plot(  [frequency frequency], [0 1.1*amplitude] , 'r--'); hold off;
            title( sprintf( 'f0=%6.3f Hz, g=%5.1f mHz, A=%5.3f mV, B=%5.3f mV',frequency, gamma*1000, amplitude*1000, base*1000) );
            xlabel('frequency (Hz)'); ylabel('Demod Amplitude (V)'); ylim([0 1.1*amplitude]);
            subplot(2, 1, 2);
            plot(phase_data.grid, phase_data.value/pi*180.0, 'b.-'); grid on; hold on;yticks(-120:30:120);
            plot(  [frequency frequency], [-120 120] , 'r--', minmax((phase_data.grid).'), [0 90 -90; 0 90 -90], 'k--'); hold off;
            legend('phase'); ylim([-120 120]); title(sprintf('slope = %5.2f deg/mHz', slope));
            
           
            
        end
    end
end

